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Abstract 

Polyharmonic spline (PHS) radial basis functions (RBFs) are used together with polynomials 
to create local RBF-finite-difFerence (RBF-FD) weights on different node layouts for spatial 
discretization of the compressible Navier-Stokes equations at low Mach number, relevant to 
atmospheric flows. Test cases are taken from the numerical weather prediction community and 
solved on bounded domains. Thus, attention is given on how to handle boundaries with the 
RBF-FD method, as well as a novel implementation for the presented approach. Comparisons 
are done on Cartesian, hexagonal, and quasi-uniformly scattered node layouts. Since RBFs are 
independent of a coordinate system (and only depend on the distance between nodes), changing 
the node layout amounts to changing one line of code. In addition, consideration and guidelines 
are given on PHS order, polynomial degree and stencil size. The main advantages of the present 
method are: 1) capturing the basic physics of the problem surprisingly well, even at very coarse 
resolutions, 2) high-order accuracy without the need of tuning a shape parameter, and 3) the 
inclusion of polynomials eliminates stagnation (saturation) errors. 


1 Introduction 

In applications of radial basis functions (RBFs) for fluid modeling (both incompressible and com¬ 
pressible), infinitely smooth RBFs have traditionally been used due to their spectral convergence 
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properties, with multiquadrics and Gaussians being the most popular [MiiiiEiiiiiinKiaiaEg. 
However, fluid flows in nature can exhibit complex rapidly developing features with such steep 
gradients that spectral accuracy can not be realized on resolutions that are observable or practical. 
This study offers a different perspective on RBF-based fluid modeling with the following aspects: 
1) using odd-ordered polyharmonic spline (PHS) RBF, for m G N (r is the Euclidean dis¬ 

tance between where the RBF is centered and where is evaluated), and 2) in conjunction with 
higher-order polynomials (degree > 4). From a historical perspective, using this combination for 
RBF-FD has not been considered, most likely for the following reasons: 

1. Before the development of RBF-FD or other flavors of local RBFs |32ll33l[37] . applications of 
RBFs were global. Then, if piecewise smooth RBFs were used, they were used in conjunction 
with low order polynomials, e.g 1, x, y in 2-D. The only role of the polynomial was to guarantee 
non-singularity of the RBF interpolation matrix, which needs to be inverted to derive the 
differentiation matrices mm- The role of capturing the physics of complicated fluid flows 
was the left to the RBFs. 

2. Even when used in a global sense, these RBFs were not as popular as infinitely smooth RBFs. 
For example, r^, results in an RBF that jumps in the third derivative, giving at best fourth- 
order accuracy in 1-D (with the order of convergence increasing as the dimension increases 
(c.f. [28]) assuming smooth data). The curse lies in the fact that as m increases, leading to a 
smoother RBF, the condition number of the interpolation matrix gravely increases. Thus, in 
the past, one was limited to keeping m small and having low algebraic accuracy. 

3. Lastly, using polynomials on a global scale can be dangerous, since it can lead to Runge 
phenomena near the boundaries. In contrast, on a local scale as in the RBF-FD method, one 
is only interested in the approximation at the center of the stencil and not at the edges. 

As a result, a new way to use PHS RBFs combined with polynomials in the context of RBF- 
FD is introduced, such that high-order accuracy is gained with excellent conditioning of RBF-FD 
interpolation matrix and no saturation error is encountered. Furthermore, there is no need to 
bother with selecting an optimal shape parameter, which plays an important role in the accuracy 
of the solution when using infinitely smooth RBFs na m El Eu Ej). We will demonstrate the 
performance of the modified RBF-FD method for 1) the advection of a scalar in a strong shear 
flow (a hyperbolic PDE introduced by [2l| and popularized by [2|) and 2) the 2D nonhydrostatic 
compressible Navier-Stokes on bounded domains applied to test cases common in the numerical 
weather prediction community. Although already broad in scope, the authors further wish to 
classify the differences, if any, that occur in applying this methodology on different node layouts: 
1) Cartesian, 2) hexagonal, and 3) scattered. The rational being that classical finite difference 
methods, based on polynomials, are usually implemented on Cartesian lattices; hexagonal layouts 
are optimal in terms of node packing in 2D, supplying information along 3 distinct directions in 
contrast to Cartesian layouts where information is aligned only along 2 directions; and scattered 
nodes allow for geometric flexibility of the domain and the ability of node refinement. 

The paper is organized as follows: Section very briefly introduces RBFs. Section discusses 
the calculation of RBF-FD weights using polynomials. Section demonstrates how the inclusion 
of polynomials with PHS eliminates stagnation (saturation) errors. Section discusses the node 
sets that are used and how boundaries and hyperviscosity are handled. Section applies the 
methodology on the various node layouts, giving detailed results from test cases that are standard 
in the numerical weather prediction community. Lastly, Section summarizes the observations of 
this paper. 
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2 A brief introduction to Radial Basis Functions 


An RBF is a d-dimensional radially symmetric function (^ : —>■ M that depends only on the 

Euclidean distance between where the RBF is centered, Xc and where it is evaluated, x. That is, 
regardless of dimension, its argument is always a scalar defined by r = ||x—Xc|| 2 , where 11-112 denotes 
the Euclidean distance. RBEs come in two flavors: 1) Piecewise smooth and 2) infinitely smooth. 
The former features a jump in some derivative and thus can only lead to algebraic convergence. The 
latter is a C°° function and can lead to spectral convergence when the data is sufficiently smooth. 
Only PHS RBEs do not depend on a parameter e that controls the shape of the RBF (which 
influences both the conditioning of the matrices and the accuracy of the results |39[ 130)1. This last 
comment is a strong motivation of this paper since, by using PHS RBF, one avoids the difficulty 
of dealing with a shape parameter and yet can achieve high-order accuracy. For theoretical aspects 
of PHS (a class of conditionally positive definite radial functions), see [111 141) . Common RBFs of 
both categories are given in Table where Xc represents where the RBF is centered. 


Table 1: Some common choices for radial functions 


Piecewise smooth RBFs 

(t>{r = X - Xc 2 ) 

Polyharmonic Splines (PHS) [9] 

r'27n g PJ 


r2m+i^ m G N° 

Matern [25] 

> 0, (Bessel A'-function) 

Compact support (‘Wendland’ flO]! 

(1 — er)'^p(£r), p certain polynomials, m € N 

Infinitely smooth RBFs 

Gaussian (GA) 


Multiquadric (MQ) 

a/ 1 -h (er)2 

Inverse Multiquadric (IMQ) 

l/\/l + /r)2 

Inverse Quadratic (IQ) 

1/(1 + (er)2) 


3 Calculation of the PHS RBF-FD differentiation weights with 
polynomials 

The differentiation weights are derived so that the resulting linear system becomes exact for all 
RBF interpolants s(x) of the form s(x) = ~ Xi|| 2 ) + te(x)} with the constraints 

— 0, where p/(x) are all polynomials up to degree I in the dimension of the problem. 
These constraints enforce that the RBF basis reproduces polynomials up to degree I as well as 
ensure that the far-held RBF expansion is regularized (i.e. does not blow up) HZ). They also are 
known as the vanishing moment conditions |23) . 

It can then be shown (see Section 5.1.4 in [18) 1 that the above leads to the following linear 
system for the differentiation weights. 


3 










I || 2 m+l 

1^1 - X 1 II 2 


I || 2 m+l 

|Xn - X1II2 


Xl 

yi 


I || 2 m+l 

1=^1 - Xn ||2 


Xr,, - X- 


12 m+l 


nil 2 


1 

Vn 


1 Xl yi 


1 Xn Vn 


0 0 0 
0 0 0 
0 0 0 


Wl 

Wn 

Wn+l 

Wn+2 

Wn+2, 


T II i| 2 m+l 

L ||x — X1II2 


L ||x — X, 


12 m+l 


n||2 


Ll 

^®lx= 

^ylx=: 


X=Xc 

X=Xc 

X=Xc 


( 1 ) 


where for simple illustration purposes, we have included only up linear polynomials. The weights 
Wn+i to Wn+ 1 , are ignored after the matrix is inverted. Solving 0 will give one row of the dif¬ 
ferentiation matrix (DM) that contains the weights for approximating L at Xc- Thus this process 
is repeated N times over all nodes in the domain, giving a preprocessing cost of 0{n^N). Since 
n « N, it should be noted that the DM usually becomes over 99% empty. As a result, we do not 
actually store the DM but only its nonzero entries. 


4 Stagnation error, PHS order, and polynomial degree 

Stagnation (saturation) error is defined as the convergence either stagnating or increasing as res¬ 
olution increases. For infinitely smooth RBFs, as the resolution increases (i.e. r decreases), the 
shape parameter must increase to maintain the condition number of the matrix in Q, resulting 
in stagnation error since the more peaked the infinitely smooth RBFs become as e increases, the 
less accurate the approximation. In contrast, PHS RBF with polynomials can achieve high-order 
algebraic convergence without encountering saturation error or the difficulty of finding an optimal 
value of the shape parameter e for good accuracy, which has been a central focus of quite a few 
studies [H [221 la n [8]. It should be noted that if polynomials are not included with the PHS 
RBF, stagnation error is encountered since boundary errors at the edge of the RBF-FD stencil can 
be quite large and penetrate toward the center of the stencil where the interpolant or any deriva¬ 
tive is being approximated. Further investigation of the effects of adding polynomials to RBF-FD 
approximations for both infinitely smooth and PHS RBFs is given in [14j . 

Since locally all smooth functions are well represented by Taylor expansions, then under refine¬ 
ment, the RBF-FD approximation must reproduce polynomial behavior. In the following numerical 
studies, it indeed was found that the convergence rate is dictated not by the order of the PHS, 
m, but by the highest degree of polynomials, I, used. In addition, for PDFs with only first-order 
spatial derivatives, as in all the test cases, the convergence rate can be expected to be 0{h^). The 
reason being is that the error for polynomial interpolation is but one order in h is lost 

in approximating a first derivative. These observations are in excellent agreement with Figure 
where d/dx of the smooth function /(x, y) = 1 + sin(4x) -|- cos(3x) -|- sin(2y) is approximated with 
two different PHS RBFs, and r^, augmented with polynomials (e.g. poly 3 is augmentation of 
the matrix in (j^ with the 10 polynomials up to degree 3 in 2D). The approximation is at the center 
of 37 node hexagonal stencil with the evaluation nodes of the derivative in its vicinity. The slopes in 
the two panels of Figure are identical, the only difference being that the constant that multiplies 
the order of convergence is slightly smaller for r^, thus giving a marginally higher accuracy for a 
given resolution. It is still important to note that the RBFs play a crucial role in safety against 
singularities due to particular node layouts. 
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Figure 1: The convergence rate in approximating d/dx of the function f{x,y) = 1 + sin(4x) + 
cos(3x) + sin(2?/) on a 37 node RBF-FD stencil based on and angmented with polynomials as 
described in the text. The dashed line marks machine round-off errors in standard double precision 
of for approximating the first derivative. 


5 Node sets, Ghost nodes, and Hyperviscosity 


Section ^ overviews the various node sets considered in this study. Section 5^ discusses how to 
increase accuracy near bonndaries with the use of ghost nodes and Section ^ discusses the need 
for and type of hyperviscosity nsed, introdncing a novel way of implementing hyperviscosity with 
PHS and polynomials. 


5.1 Node-sets 

Unlike traditional methods, RBF-FD has the advantage of being equally simple to apply on any 
set of nodes. Figure shows the three types of node-layouts that will be considered in the present 
study. First, we consider Cartesian since they are the lattices classical finite differences (FD) are 
usually implemented on. Next, it is well known that hexagonal node sets are the optimal packing 
strategy in 2D (i.e. for a fixed area, one can ht the most number of nodes). They have also 
been considered an optimal layout for differentiation stencils, since information is aligned along 
three different directions as opposed to only two with Cartesian layouts. Although FD have been 
sporadically implemented over the decades on hexagonal node layouts, they have never caught 
favor do to the complexity of the implementation. However, now with RBF-FD, implementation is 
simple. Lastly, scattered node layouts are considered as they have the great advantage of geometric 
flexibility that will be needed in fntnre applications of RBF-FD on irregular domains and/or with 
local node refinement. 
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Scattered 


Figure 2: Three different types of node-distributions that are used to solve the 2D test problems. 

5.2 Ghost nodes 

Near boundaries, stencils become one-sided, leading to a deterioration of the approximation due 
to Runge phenomenon. In order to ameliorate this effect, one layer of ghost nodes can be placed 
just outside the domain. Once the ghost nodes are placed, function values at these locations can 
be solved for by enforcing additional constraints at the boundary nodes. For example, if the upper 
boundary of a rectangular domain is free-slip, then dujdz = 0, and this condition can be used to 
solve for values of u at the ghost nodes. Similarly, the PDF itself can be enforced on the boundary, 
giving an extra constraint. 


Interior Node 

• 

Boundary Node 

• 

Ghost Node 

o 

Evaluation Point 

• 

Stencil Boundary 
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Figure 3: An example of an RBF-FD stencil that might be used for enforcing dujdz = 0 on the 
top boundary. 


For each node on the top boundary, one ghost node is placed just outside the boundary, as shown 
in Figure]^ The function values for u at the ghost nodes are obtained by enforcing dujdz = 0 
at each of the top boundary nodes simultaneously. This will lead to a coupled system with as 
many equations as there are ghost nodes. In the case illustrated in Figure the system will be 
tridiagonal, as there are three unknown values at the ghost points for each evaluation node on the 
boundary. 

The following is a more detailed discussion on how the ghost node values are calculated. Suppose 
there are Nj interior nodes, Nb boundary nodes, and Nq ghost nodes, and let the total number 
of nodes be = Ni Nb + Nq- For each top boundary node, approximate the differentiation 
weights for djdz as given in Section]^ This will result in a sparse Nb x N DM (here, called W). 
Thus, 

du ^ , 

=0 IS approximated by 
dz 

Wu =0. 
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The function values u are organized according to where they are located: 

Uj : function values at interior nodes 
: function values at boundary nodes 
Uq : function values at ghost nodes 

Then, the condition Wu ~ 0 can be written as 

lILi T ~ Oj (^) 

where the matrix W has similarly been split into pieces according to the three different types of 
nodes: 


i^B X : weights applied on interior nodes 
Wb {Nb X Nb) : weights applied on boundary nodes 
Wg {Nb X Ng) : weights applied on ghost nodes 

Finally, (§ is used to solve for Uq, the function values at the ghost nodes: 

Uq = —bTg ^ {W^jUj + WbUb) • 

Once the function values at the ghost node are known, they can be used for the approximation of 
other derivatives that appear in the governing equations. 

5.3 Hyperviscosity with PHS and polynomials 

When the viscosity of the fluid /x is small (such as the case with air ~ 10“®m^/s), there is essentially 
no natural diffusion in the governing equations, and high-frequency errors will grow to dominate 
a numerical solution. To achieve time stability with the RBF-FD method, it has been shown that 
adding a relatively small amount of hyperviscosity to the right-hand-side of the governing equations 
eliminates the contaminating high-frequency noise while keeping the numerically relevant portion 
of the solution intact [HKiais]. 

The hyperviscosity operator takes the form yA^, where k is the power the Laplacian and 7 is 
a scaling parameter. The integer k controls which frequencies are most affected, with larger values 
of k giving stronger damping of high frequencies and weaker damping of low frequencies. As has 
been shown in US das], for good stability and accuracy, the parameter 7 is directly proportional 
to the number of nodes in the domain N or conversely the resolution /i, as well as k. Thus, for a 
square-type domain in 2D h ~ 1/y/N, and 7 = where c is a constant that is generally set 

for the problem at hand and independent of the resolution h, (e.g. for the NS test cases c = 
regardless of the resolution or node layout used). 

The hyperviscosity operator is particularly simple to apply to an odd-powered PHS RBF, 
regardless of the spatial dimension. The Laplace operator in d dimensions for a radially symmetric 
function is given by A = jdr^ -|- {{d — l)/r)d/dr. Apply this to 4>{r) = r™ results in 

A (||x|r) = m[m + {d- 2)] Ijxir'^ . (3) 

In other words, applying the Laplace operator to a PHS RBF of degree m gives a new PHS RBF of 
degree m — 2. Using the above relationship, one can evaluate A^ (||x||™'), and the new RBF will be 
continuous provided that (m > 2A: -|- 1). Higher-order Laplacians can be implemented by applying 
Q repeatedly. It should noted that the order of the PHS used for spatial discretization need not 
be that used for hyperviscosity. However, it was found experimentally that the simplest approach 
for the needed inclusion of polynomials was to use up to the same degree as that for discretization. 
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6 Numerical Studies 


The first test case, inviscid transport of a scalar variable in a strongly sheared vortex flow, is a 
case of pure advection with a known analytical solution, so that the convergence properties of the 
method can be tested. It was originally proposed in |24j and then considered in the context of 
applying limiters in 121135]. The second set of tests is based on the work presented in [38], where 
a cold descending bubble in a neutrally-stratified atmosphere develops into a traveling density 
current with the formation of Kelvin-Helmholtz rotors. It is now considered a classic test case in 
nonhydrostatic atmospheric modeling. The third test [29] (with similar studies in |2I]) simulates a 
rising thermal air bubble and nicely illustrates how instability patterns at the leading edge of the 
thermal are dependent on the node layout when the the dynamic viscosity is that of air. 

In all cases, n = 37 node stencils are used for spatial discretization, since both Cartesian 
and hexagonal layouts have perfectly symmetric stencils at that number, as seen in Appendix A. 
Although not essential, stencil symmetry is beneficial in that it it provides information evenly for 
approximating an operator at the center of the stencil. In all cases, time stepping is done with a 
4th-order Runge-Kutta scheme (RK4). 


6.1 Advective transport of a scalar variable 


In this test, a circular scalar field is stretched and deformed into a crescent by vortex-like velocity 
field that then reverses and returns it back to its original position and shape. The governing 
equation is defined on [0,1] x [0,1] in x and y and given by 


dip 




A 

dy 


{vTp) . 


The scalar pj is advected by the following divergence-free velocity field 


u (x, y, t) = uq (r, t) sin 6, v {x, y, t) = —uq (r, t) cos 9. 

with period T, where 


dvrr 

U0 (r, t) = ^ 


1 — cos 


27rA 1 — (4r)^ 


and 


T J 1 + (4r)6j ’ 
r = \J{x- 0.5)^ + (y - 0.5)^ 9 = tan“^ ' 


In order to create a test problem with no boundary effects, the nodes on the interior of the 
domain near the boundary, as well as the function values associated with them, are simply reflected 
over the boundary (thus no boundary conditions are needed), forming perfectly symmetric boundary 
stencils - half of which are then ghost nodes and half interior nodes. The boundary is then time 
stepped with the rest of the interior of the domain. An example of a symmetric 19 node boundary 
stencil is given in Figure]^ 

The initial condition for ^ is a cosine bell 


V’li=0 = 


l+cos(7rr) 

2 

0 


r < 1 
r > 1, 


(4) 


where r = Sy (x — 0.3)^ + {y — 0.5)^. Using a classic Runge-Kutta 4th-order (RK4), ip is advanced 
in time from t = 0 until t = T (one period), at which point it should have ideally returned to its 
original height and position. 










Figure 4: An illustration of a hexagonal 19 node boundary stencil in which the nodes and there 
associated function values are reflected outside the domain. This setup is used for the advection of 
a scalar transport test case. 

6.1.1 Accuracy of solution, convergence, and effect of different node layouts 

Figure shows the time evolution for solution (assuming a period of T = 1), with corresponding 
velocity field, on a hexagonal node layout. The specifications of the resolution, time-step, basis 
functions and hyperviscosity used are given in the caption of the hgure. Although, this is a high 
resolution case with a total of A = 40401 nodes in the domain or a node spacing of h = 0.005, 
it should be noted that the maximum amplitude of V’ has only increased by 0.07% and has only 
gone below zero by —0.08%, as can be seen in Figure In fact, when the scalar held is in its 
state of highest deformation at t = 0.25, 0.75, the error in the maximum amplitude is not more 
than 0.001. Furthermore, the CFL criterion for an RK4 stability domain dictates a time step of 
At < Ax • 2\/2/(maxu), which for this case translates into At < 7.0(10)“^, a factor of only 2.12 
larger than the time step taken of At = h/15 = 0.005/15 = 3.3(10)“^. 

In order to observe the long-time errors in the method, the solution is advanced for 100 periods 
as shown in the left panel of Figure Even after so many revolutions, the height of the scalar 
held has decreased only by 4%, with a slight distortion from its circular shape. The right panel of 
Figure shows a dispersive error pattern with a maximum value of 0.13. The £2 and £00 errors are 
0.125 and 0.141, respectively. 

For a 37 node stencil, the highest degree polynomials that can be used on all three node sets is 
hfth degree. Both hexagonal and scattered node sets can handle sixth degree polynomials but not 
Cartesian layouts. This is because on such a lattice the nodes approach non-unisolvency, resulting 
in the column vectors of the polynomial portion of the matrix in Q becoming linearly dependent. 
Also, in order to demonstrate that the polynomial degree controls the convergence and not the 
PHS order, PHS are now used instead of r® with up to hfth-degree polynomials on Cartesian, 
hexagonal and scattered nodes. Figure [^illustrates this for three resolutions, h = 0.02,0.01,0.005. 
For any given resolution, all nodes sets perform roughly the same, both with regard to the minimum 
and maximum function values and the errors in the £2 and £00 norm. Comparing the the maximum 
and minimum of the solution for hexagonal nodes on h = 0.005 between Figures andit can been 
seen that using r® gives slightly better accuracy. This phenomena was noted in Section [^ Figure 
shows the convergence rate corresponding to the cases given in Figure From the discussion in 
Sectionj^ the convergence rate should be 0{h})^ where I is the highest degree of polynomials used. 
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Figure 5: Time series of the solution for if; at: (a) t = 0(1), (b) t = 0.25, (c) t = 0.5, (d) t = 0.75, 
with the corresponding minimum and maximum values at each time. Contour lines are in intervals 
of 0.05. An n = 37 node stencil with r® PHS and up to 4*^-order polynomials on a hexagonal 
node layout of 201 by 201 is used. The time-step is h/15, where here h = 1/200 = 0.005. A 
hyperviscosity of —is also implemented. 


in this case 5. Fifth-order convergence is indeed seen for all node sets in Figure Also in this 
figure, the RBF-FD method is compared to a 5th-order upwind scheme with and without a WENO 
limiter, both of the latter exhibiting a third-order convergence rate. The reason for the comparison 
is that this order FD upwinding scheme is the type used in the Weather Research and Forecasting 
(WRF) Model (http://www.wrf-model.org/). The time steps for both methods are comparable, 
with less than a 1% difference. 

It should be remembered that this test case was set up to investigate what the numerical results 
for the proposed RBF-FD method would be under no boundary effects. So, if there is the unusual 
circumstance that the solution does not interact with any boundaries in a bounded domain and no 
refinement will be needed, then solving the problem on a Cartesian lattice will give just as good 
results as hexagonal. Furthermore, the fact that scattered nodes performed just as well as the 
other two layouts is of great benefit since it paves the way for the ability to implement local node 
refinement. 
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Figure 6: Solution (left panel) and error (right panel) at t = lOOT on 40,401 {h = 0.005) hexag¬ 
onal nodes using (j){r) = r® with up to 4*'^ degree polynomials on a 37-node stencil and A^-type 
hyperviscosity. 
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Figure 7: Plots of solutions at t = T = 1 using the three different types of node-sets. </>(r) = 
polynomials up to degree 5 on a 37-node stencil and A^-type hyperviscosity were used. The amount 
of hyperviscosity, 7 , varies between node-sets, but is on the 0 ( 10 )“^^ to 0 ( 10 )“^“^. 


12 







































Figure 8; Convergence plots of relative error ^^epctlla ^ (/i = 0.02,0.01,0.005, 0.0025). 

In all cases, (f){r) = with polynomials up to degree 5 on a 37-node stencil and A^-type hypervis¬ 
cosity were used. The error decreases O (/i^), which is expected since up to 5**^ degree polynomials 
were included and only first derivatives need to be approximated in the PDF. 

6.2 Governing Equations for Navier-Stokes test cases 

In all test cases below, the set of governing equations is the 2D nonhydrostatic compressible Navier- 
Stokes equations at low Mach number, M~ 0.1, in a rectangular or square domain. The equations 
are given by 


du 

du 

du 

m ^ 


~'"Yz ■ 

dw 

dw 

dw 




dd 

d6 

de 

m ~ 



dir 

dix 

diT 

~di ^ 




a— 

a— 


Rd 


-TT 


- g + gAw, 


du dw\ 
dx^ dz ) ' 


where u and 

p \ RdICp 


PoJ 


w are the velocities in the horizontal and vertical directions, respectively, tt = 
is the non-dimensional Exner pressure (Pq = 1 x 10® Pa); and 0 = ^ is the poten¬ 
tial temperature. The constants Cp = 1004 and Cp = 717 are the specific heat at constant pressure 
and the specific heat at constant volume, respectively, with the gas constant for dry air being 
Rd = Cp — Cp = 287. Additional parameters are g = 9.81m/s^ , the gravitational constant, and g,, 
the dynamic viscosity. Furthermore, it is assumed that all quantities to be solved for, [u,w,6,ttY , 
are perturbations {') to a background state (“) that is in hydrostatic balance, i.e. the fluid is initially 

at rest, u = w = 0, and the background Exner pressure is a linear function of height z, — = — 

dz CpO 

Substituting this latter relation into the equations above and writing 9 = 6 + 6' and tt = n + tt' 
(where the (') symbol has been dropped for reading clarity) yields the governing equations to be 
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used for computation: 



du 

du 

du 

m 



dw 

dw 

dw 


~'^~di 


do 

de 

de 

m 


— w ——h 
dz 

Ott 

d-K 

( dw 





Ott 
dx 
. dTT 

’ Ih 


gO 

6 


+ 


dz 




(5) 

( 6 ) 

(7) 

( 8 ) 


In the following studies, the perturbation notation (') is generally included when reporting results 
to keep in mind these are perturbation quantities. 


6.3 Numerical Set-up for the NS cases 

The governing equations ([^-Q are solved numerically using a method-of-lines (MOL) approach. 
PHS RBFs, (f){r) = with polynomials up to fourth degree are used on a stencil-size of n = 37 
to approximate all spatial derivatives locally. The remaining system of first order ODEs in time is 
solved with RK4. A A^-type hyperviscosity is applied in all cases to damp high-frequency modes. 
The time-step for all test cases as a function of resolution is 


node-spacing (h): 

800m 

400m 

200 m 

100 m 

50m 

25m 

At: 

2 s 

Is 

is 

2 ® 



As 

16® 


Table below gives the domain size, and number of nodes used as a function of resolution for 
the numerical studies of the NS test cases. 


Table 2: Information regarding the computational domain for each test case. The number of nodes 
{N) is for hexagonal nodes. 


Test (domain size (x x z in km)) 

Resolution (m) 

« N 

Straka Density Current [38] ([—25.6,25.6] x [0,6.4]) 

800 

720 


400 

2,700 


200 

10,000 


100 

38,500 


50 

152,650 

Translating Density Current [42] ([0,36] x [0,6.4]) 

800 

500 


400 

1,900 


200 

7,040 


100 

27,040 


50 

107,350 

Rising Thermal Bubble ([0,10] x [0,10]) 

200 

2,980 


100 

11,760 


50 

46,720 


25 

185,430 
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6.4 Density Current 

In the density current test case [38], a hydrostatic neutral atmosphere is perturbed by a bubble 
in the potential temperature. A mass of cold air falls to the ground and develops three smooth and 
distinct Kelvin-Helmoltz rotors as it spreads to the sideways. This test has become widely used 
in weather modeling community for assessing the ability in new numerical schemes to capture the 
physics in nonhydrostatic fluid flows |27[ 1201 [36lI26| . Figure]^ shows the behavior of the numerical 
solution in time from t = Os until the final time, t = 900s. 

The computational domain is [—25.6,25.6] x [0,6.4] km^, and the governing equations Q-Q 
are solved with a viscosity of /r = 75 m^/s. 

Define 9 and If 

Let T{z) = Ts — -^zhe the background state for temperature, where Tg = 300 is the temperature 
at the ground surface in Kelvin. Then, the background states for potential temperature and Exner 
pressure are given by 


^ ^ 1 9 

= Ts, Tr{z) = - =1-— 

^ ^ 9 c„Tg 


z. 


Define initial conditions 

The vector of unknowns is initially zero except for the potential temperature. 

u\t=o = 0> w\t=o = ^'L=o = 0- 

The (C^) initial condition for 9' is derived via a cool cosine bubble in the temperature T defined 
by 


r'l = 

h=o 


— ^ {1 + cos [7rcr(x, 2;)]} , r{x,z)<l, 


0 , r{x, z) > 1, 
where tTc = 3.14159... is the standard trigonometric constant and 


(9) 


r{x, z) = 


X — X, 


Xr 


+ 


Z — Zr 


(xc, Zc) = (0 km, 3 km), 
(xr, Zr') = (4 km, 2 km). 


Then, the initial condition for 9' can be found by dividing by vr: 

«'Uo 


T 

-Tg = 

T + T' 

-Tg = 

T'l 
9 + — 

TT 

t=o 

vr 

t=o 

TT 


-T, = 


it=o 


J t=o 


TT 


Define bonndary conditions 

The problem is periodic in the x direction with the following conditions on the top and bottom 
boundaries in z: 


w = 0 , 




These are the only boundary conditions necessary to solve the governing equations. However, 
enforcing the vertical momentum equation Q on the top and bottom boundaries and assuming 
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that perturbation in the pressure gradient balances the perturbation in the potential temperature 
leads to the following condition for tt' , 

dn' _ g9' 

~ CpO {9 + 9') ■ 

Furthermore, since the dynamic viscosity g for air is ~ 10“^, a good approximation on the top and 
bottom boundaries is Aw' = 0 or d‘^w'(dz^ = 0 since rc = 0 on these boundaries. While these two 
extra boundary conditions on tt' and w' conditions are not required, they allow for the use of ghost 
nodes in all four variables. In summary, the lateral boundaries are periodic, and the complete set 
of boundary conditions enforced on the top and bottom boundaries is; 

_d^w _du _ d9' _ d^' _ g9' 

^ ~ ~ dz ~ Ih ~ ’ Ih ~ Cp9 {9 + 9')' 



Figure 9: Time evolution of the potential temperature perturbation 9' for the density current test 
case with g = Ihrr? js at a 100m resolution on hexagonal nodes. 
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Figure 10: Convergence behavior for O' in the Straka density current test case with /x = Ibrn^/s. 
The h =800m, 400m, 200m, 100m, and 50m errors were calculated using the 25m RBF-FD reference 
solution. 


Cartesian Hexagonal Scattered 



-12.5-10.5 -8.5 -6.5 -4.5 -2.5 -0.5 -12.5-10.5 -8.5 -6.5 -4.5 -2.5 -0.5 -12.5-10.5 -8.5 -6.5 -4.5 -2.5 -0.5 


Figure 11: The solution at the final time for the density current test case solved on different node 
layouts for resolutions 800m, 400m, 200m, 100m. Only half the solution is shown to enlarge details. 
Contours for the density current begin at -0.5K and are in intervals of IK. The white areas are 
enclosed by a contour of 0.5K. 
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As can be seen in Figure 10 the error in the £2 norm is approximately the same regardless of 
the node layout with a third-order convergence rate for node resolutions 200m and less. However, 
convergence rates do not illustrate how the physics is being resolved with regard to where the data 
is sampled (i.e. in terms of the node layout). As a result, in Figure [TH the solution of the density 
current test case is given on the three different node layouts discussed in Section for four different 
resolutions, varying from 800m to 100m. In the highest resolution displays (100m), although all 
node layouts seemed to have converged to the same solution, differences can be noted in Table 
where a 25m test run on hexagonal nodes is used as a reference solution. Results from a 50m test 
run on the 3 different node sets are also given in the table. On hexagonal and scattered nodes the 
minimum 0' has indeed converged by 100m, while the maximum 9' is still a fifth of a degree off. 
Similar error percentages can be found in w' and in the front location at these fine resolutions of 
100m and 50m, noting that Cartesian nodes perform slightly worse. Notice that the front for the 
100m Cartesian is at the same location as achieved by a 200m hexagonal layout. 

At coarser resolutions, such as 400m and 800m, values in the table will be far off the converged 
25m solution. Instead, noting physical features of the solution in Figure 11, such as 1) at what 
resolution do the rotors begin to form, their shape and where, and 2) how much cold air has been 
entrenched in each rotor, will give a better idea of the capability of the node layout to capture the 
physics. The following observations can be made: 


1. At 800m - approximately 720 nodes in the domain: The hexagonal and scattered node calcu¬ 
lations give more clear evidence of the first (largest) rotor being formed. The -3.5K contour 
in the hexagonal case (circular inner most contour in the first rotor) is even close to its final 
position if compared to the 100m case. Although the first rotor for the scattered case is not 
quite as nicely formed as in the hexagonal case, it has entrenched more cold air, having a 
-4.5K contour (teardrop shape). Notice that at 100m the -4.5K contour is the coldest that 
appears. In comparison, the 800m Cartesian has barely any rotor formation and is much 
more wildly oscillatory, which is also noted by the fact that the maximum O' is 2.43K, at least 
1.3K larger than for the other node layouts. See Table Also note in the table that the error 
in the front location decreases from 4% to 2%, when hexagonal nodes are used opposed to 
Cartesian, with scattered given an intermediate error of 3%. Both hexagonal and scattered 
nodes undershoot the correct position while Cartesian overshoot it. 

2 . At 400m - approximately 2700 nodes in the domain: At this resolution, oscillations due to 
boundary error effects especially in the first rotor are very pronounced on the Cartesian layout; 
this carries over even to the 200m resolution for this node case. For scattered nodes there are 
minor oscillations in the solution. For the hexagonal case, barely any are evident. Formation 
of the second rotor has the nicest intact shape with the least amount of oscillation in both 
the hexagonal and scattered, with the latter having entrenched a slightly larger amount of 
cold air (notice the size of the -3.5K contour teardrop-shaped area in the second rotor of the 
scattered case). 


The differences between the columns of subplots reflect only the intrinsic resolution capabilities 
of the different node layouts for capturing the physics. The traditional Cartesian choice is the least 
effective one. At every resolution level, the hexagonal and scattered choices give better accuracy 
than the Cartesian one. The advantage of generalizing from hexagonal to quasi-uniformly scattered 
nodes, is that it then becomes easy to implement spatially variable node densities, i.e. to do local 
refinement in select critical areas. It is very important to note that this major increase in geometric 
flexibility (from hexagonal to quasi-uniformly scattered) hardly has any negative effect at all on 
the accuracy that is achieved, nor on the algorithmic complexity of the code. 
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Table 3: Resolution (h), minimum and maximum values for 9' and w, and front location at various 
resolutions for the density current test case with /r = 75m? js. The front location was determined 


.5K contour 

ine. 
h (m) 

min {9'} 

max{0'} 

min{r(;'} 

max{r(;'} 

front (m) 

Cartesian 

800 

-7.74 

2.43 

-9.19 

11.00 

16,079 


400 

-13.45 

1.10 

-15.21 

16.36 

16,013 


200 

-12.15 

0.57 

-16.59 

17.49 

15,799 


100 

-9.84 

0.27 

-16.14 

13.45 

15,500 


50 

-9.71 

0.04 

-15.96 

12.86 

15,424 

Scattered 

800 

-8.60 

1.11 

-10.11 

10.05 

15,477 


400 

-12.03 

1.13 

-13.26 

12.79 

15,747 


200 

-10.40 

0.42 

-15.90 

13.60 

15,597 


100 

-9.70 

0.21 

-16.00 

13.12 

15,447 


50 

-9.70 

0.02 

-15.95 

12.87 

15,422 

Hexagonal 

800 

-6.90 

1.00 

-11.53 

9.61 

15,101 


400 

-13.38 

0.98 

-12.93 

10.11 

15,721 


200 

-11.42 

0.44 

-15.90 

14.34 

15,501 


100 

-9.70 

0.20 

-15.90 

12.96 

15,444 


50 

-9.70 

0.01 

-15.93 

12.90 

15,420 

Reference 

25 

-9.70 

0.00 

-15.93 

12.90 

15,418 


6.4.1 Low-Viscosity Density Current /r = 2 x 10 ^ m^/s 


Here, the density current test case is repeated, except with the dynamic viscosity /r set to that of 
air. The purpose of this test case is to show that one can stably time step the RBF-FD method in 
a completely turbulent regime. The same amount of hyperviscosity as well as the same time step 
are used in this test case as in the one with n = 75 m? js. Time stability is governed solely by the 
fact that the time step could not exceed the speed of sound in air. 

At such low viscosity, the solution enters the turbulent regime. In such regimes, there is no 
convergence to any solution as energy cascades to smaller and smaller scales, eventually entering 
the sub-grid scale domain. Nevertheless, it is interesting to observe whether the model remains 
stable in this regime. Figure shows the solution at 100m, 50m and 25m resolutions on the three 
different node layouts. For any given resolution the solution looks completely different depending 
on the node layout. This is to be expected as changing the node layout in practically the absence 
of explicit viscosity is equivalent to introducing slight perturbations in the solution. A more robust 


illustration of this will be given in the test case of a rising thermal bubble. Section 6.5 
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Scattered 



Figure 12: Potential temperature perturbation O' for the low-viscosity density current 6.4.1 using 
100m, 50m, and 25m nodes at t = 900s. 


6.4.2 Translating Density Cnrrent, /r = 7hm?js 


This test is the same as in Section 6.4, except that the domain is now [0, 36] x [0, 6.4] km^ and there 


is a horizontal background wind of u = 20 m/s. The size of the domain is set up so that at t = 900s 
the two “halves” of the solution should be symmetric about x = 18km. The introduction of a 
background mean flow introduces a large difference in the movement of each half of the solution. 
The right portion of the outflow has horizontal velocities ~ 50 m/s, while the left portion has 
velocities ~ 10 m/s. As a result, it tests the ability of the scheme to translate the features of the 
solution at the correct speeds and to generate the correct rotor structures that arise from the local 
shearing instabilities. For this case, only Cartesian and hexagonal nodes are considered, as the case 
tests the degree to which symmetry is broken in the two halves of the solution at the final time. 
Figure 13 illustrates the time series of the potential temperature 9' held, showing how the right 
part of the solution is advected through the the right side of the domain, with the front locations 
facing one another at 900s (instead of facing the lateral boundaries as in the previous test case). 

In Figure 14 the two halves of the solution are compared about the line of symmetry (18km) 
for Cartesian and hexagonal nodes from 800m to 100m resolutions. The general observations that 
can be seen are: 


1. The 800m hexagonal node layout performs highly superior to the Cartesian both in terms of 
symmetry between the two sides, intactness of the large rotor, and its relative location when 
compared to the 100m case. 

2 . At both the 400m and 200m Cartesian case, the left half of the solution that has been advected 
through the right boundary displays a significant amount of Runge phenomena (‘wiggles’ 
in the contour lines near the boundary at 200m and severe distortion of the primary and 
secondary rotor at 400m). This is not the case for hexagonal nodes, which at 400m and 
200 m, shows relatively nice symmetry between the two sides. 
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3. At 100m, there is no distinction between the two node sets. 

The actual front locations, in terms of there distance from the 18km mark, are given in Table 
Note that the front on the right is farther from the line of symmetry (18km mark) than the one 
that has been advected through the boundary for the resolutions 800m to 100m. At 50m, both the 
Cartesian and hexagonal case is symmetric about the 18km mark. However the distance from the 
front to the line of symmetry varies between the two cases, 2586m versus 2595m, respectively. 

-16.5 -15.5 -14.5 -13.5 -12.5 -11,5 -10.5 -9,5 -8.5 -7.5 -6.5 -5,5 -4.5 -3.5 -2,5 -1.5 -0.5 



°0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3 3.3 3.6 


Figure 13: Time evolution of the potential temperature 9' for the translating density current test 
case. Snapshots were generated using the 100m RBF-FD solution on hexagonal nodes. 
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Table 4: Left and right front locations for the translating density current test case 6.4.2 as given 


by the distance from the —0.5K contour line to the center of the domain, x = 18km. All values 
are calculated at the hnal time, t = 900s. Results are for (/>(r) = with up to fourth degree 
polynomials on a 37-node stencil. A 25m reference solution is given for the hexagonal nodes. 



h 

left front (m) 

right front (m) 

Cartesian 

800 

3,065 

3,215 


400 

1,915 

2,085 


200 

2,098 

2,205 


100 

2,487 

2,512 


50 

2,586 

2,586 

Hexagonal 

800 

3,410 

3,575 


400 

1,956 

2,046 


200 

2,165 

2,260 


100 

2,555 

2,580 


50 

2,595 

2,595 

Reference 

25 

2,595 

2,595 


6.5 Rising Thermal Bubble 

With this last test case, the paper comes full circle in that the presented RBF-FD method is 
tested on a problem with sharp gradients and very little boundary interaction as in the advective 
transport of a scalar variable, but is modeled by the same 2D nonhydrostatic compressible Navier- 
Stokes equations as in the density current tests. The only difference is the initial condition is given 
by a cone-shaped perturbation. The bubble is warmer than the surrounding atmosphere and 
thus rises toward the top boundary. However, the time interval and domain size are chosen so that 
the bubble never interacts with the boundaries. 

There are two variations for this test case: 

1. ^ = 10 m^/s: By adding a small amount of explicit viscosity, the convergence behavior of the 
solution can be studied. 

2. ^ = 2 X 10“® m^/s: At such low viscosity, the bubble is in a turbulent regime, and the behavior 
of the true solution is unknown. The purpose of the test is to see if the RBF-FD method can 
give reasonable results when the initial condition is not even continuously differentiable as 
well as observe how the instability pattern at the leading edge of thermal changes with the 
node layout. 

6.5.1 Case /r = 10 m^/s 

The computational domain is [0,10] x [0,10] km^. The hydrostatic background states are defined 
hy 6 = Ts and ^{z) = 1 — y^z:, with Tg = 300K being the surface temperature. The horizontal and 
vertical velocities and the Exner pressure perturbation (vr') are initially zero, while the potential 
temperature perturbation is prescribed as a warm cone-shaped “bubble” with a jump in the first 
derivative (C^): 

= 2max{0,1 — r(x, z)/R} . 
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Here, R = 1.5 km is the radius of the bubble, and 


r(x, z) = \J{x - Xcf + {z- Zcf, (xc^c) = (5 km, 3 km). 

The same boundary conditions as in the density current problem are enforced on the top and 
bottom boundaries 


92 


w 


w = 




du 

dz 


de^ 

dz 


= 0 , 


The lateral boundary conditions are given by 


u 


u = 


dx^ 


dw 

dx 


dx 


di:' 

dz 


d7r' 


9O' 


cpO {e + O') ■ 


( 10 ) 


dx 


Figure [I^ shows the time series of the solution for a 25m resolution {N = 185, 730) on hexagonal 
nodes using with up to 4th-order polynomials and a A^-type hyperviscosity. The main purpose 
of this test is to make sure that, under refinement, all node layouts converge to the same solution 
(as this will not be the case in the next variation of the test) and to see if the convergence rate 
follows the predictions of Section even with a initial condition. Figure 16 shows the final 
solution for the three different node layouts from a resolution of 200m to 25m. At 25m resolution, 
all solutions are visually identical. At coarser resolutions, as 100m, Cartesian and hexagonal nodes 
are more similar with scattered nodes having more incongruities at the leading edge of the rising 
bubble. A possible reason for this could be symmetry-breaking associated with scattered node 
layouts that would affect areas of large shear. 

In terms of convergence, Figure [T tI shows that even though the initial condition is C*^, the method 
does achieve 4th order convergence as predicted when using up to fourth degree polynomials. All 
nodes sets converge at fourth order under refinement, with Cartesian giving the best accuracy for 
h < 100m. Note that at coarser resolutions from 400m too 100m only second order convergence is 
achieved. 

Table shows to what degree the solution has converged in terms of how high the bubble should 
have risen, 0', and w'. Note that when comparing against the 12.5m hexagonal node reference 
solution, all solutions in all variables have converged by 25m. In terms of max{0'} and maxjrc'} all 
solutions have converged by 50m. The min{0'} is almost identical for all node sets while min{r(;'} 
varies between node sets for coarser resolutions. In terms of bubble height Cartesian nodes seem 
to perform the best. 



Figure 15: Time evolution of the potential temperature 9' for the /r = 10 w?/s rising thermal 
bubble. Snapshots were generated using the 25m RBF-FD solution on hexagonal nodes. 
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Figure 16: Numerical solutions for the rising thermal bubble with /i = 10 m^/s (6.5.1) on the three 
different types of node distributions at various resolutions, shown at the final simulation time, 
t = 1100s. 
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Figure 17: Convergence behavior for O' in the rising thermal bubble test case 6.5.1 The h =400m, 
200m, 100m, 50m, and 25m errors were calculated using the 12.5m RBF-FD reference solution. 


Table 5: Resolution (/i), minimum and maximum values for O' and w', and bubble height at various 
resolutions for the rising thermal bubble 6.5.1 Results are for (j){r) = with up to fourth degree 


polynomials on a 37-node stencil. The bubble height was determined by the intersection of the 


O.IK contour and the line x = 5km. 



h (m) 

min {O'} 

max {O'} 

minjtc^} 

maxjtc'} 

bubble height (m) 

Cartesian 

200 

-0.11 

1.46 

-7.56 

11.06 

8,467 


100 

-0.08 

1.53 

-7.93 

11.34 

8,539 


50 

-0.02 

1.43 

-7.87 

11.43 

8,534 


25 

0.00 

1.43 

-7.74 

11.43 

8,535 

Hexagonal 

200 

-0.11 

1.36 

-7.67 

11.12 

8,686 


100 

-0.09 

1.65 

-8.05 

11.49 

8,527 


50 

-0.02 

1.43 

-7.75 

11.43 

8,553 


25 

0.00 

1.43 

-7.74 

11.43 

8,535 

Scattered 

200 

-0.11 

1.25 

-7.74 

10.92 

8,581 


100 

-0.09 

1.48 

-8.42 

11.40 

8,557 


50 

-0.02 

1.43 

-8.22 

11.43 

8,525 


25 

0.00 

1.43 

-7.75 

11.43 

8,535 

Reference 

12.5 

0.00 

1.43 

-7.74 

11.43 

8,535 


6 . 5.2 Case /r = 2 x 10"^ m^/s 

The rising thermal bubble test case is repeated with the viscosity n set to 2x 10“^. The first purpose 
of this test case is simply to demonstrate that the proposed RBF-FD method, implemented with 
such low viscosity and a initial condition, has complete time stability using the same time step 
and amount of hyperviscosity as is the previous section. Secondly, we are interested in observing how 
the instability pattern at the leading edge of thermal bubble evolves as the node layout changes. 
Normally, in numerical testing, to see different evolutions of a solution the initial condition is 
perturbed. However with RBF-FD, one has the flexibility of leaving the initial condition intact 
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and perturbing the node layout, which in a turbulent regime will lead to different evolutions of the 
solution. This can be seen in Figure [Ts} At 200m, there is not much difference between the bubbles. 
However as can be seen in the 25m results, the shear instability layer at the leading edge of the 
bubble (darkest contours) develops tight eddies whose structure varies significantly depending on 
the node layout. In both the Cartesian and hexagonal case, the eddy development is completely 
symmetric about the midpoint of the bubble due to the symmetry in the node layout, while in the 
scattered node layout this is not the case (a seemingly more realistic scenario for modeling warm air 
entrainment in the atmosphere). Furthermore, the scale of the eddies and the degree to which they 
excite finer scale instabilities varies between the node sets. The Cartesian node layout produces 
the largest scale eddies as well as the smallest amount of eddies. In contrast, the hexagonal nodes 
produce a rather strange bubble shape with very fine scale eddy structure. 
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Figure 18: Numerical solutions for the rising thermal bubble with fj, = 2 x 10“® m^/s on the 
three different types of node distributions at various resolutions. All results are shown at the final 
simulation time, t = 1100s. 

7 Conclusions and summary 

In this paper, a modified RBF-FD method is introduced to construct differentiation weights based 
on a combined RBF-polynomial basis, using RBF polyharmonic splines {(p{r) = r^) with polyno- 
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mial functions up to degree I in the given dimension of the problem. The method is applied to 
three standard test cases in the numerical weather prediction community [21 [38l 129] , with the latter 
two based on the Navier-Stokes (NS) equations. In addition, the effect of node layout (Cartesian, 
hexagonal, or scattered) on the error as well the qualitative character of the solution is considered. 
The following observations are made: 

1. Under refinement, with the inclusion of polynomials, stagnation (saturation) error is evaded. 

2. In the absence of boundary effects, the convergence rate is controlled, not by the order of the 
PHS, but by the highest degree polynomials used. 

3. Increasing the order of the PHS marginally increases the accuracy, as the constant that 
multiples the convergence rate decreases. 

4. For stable configurations that require no tuning of the hyperviscosity (e.g. 7 = for 

all NS tests), the number of nodes in the stencil, n, should be approximately twice the number 
of polynomial basis functions, (/+l)(/+2)/2 in 2D. Hence on an n = 37 node stencil, up to 
fourth-order polynomials (15 in 2D) are used for the NS equations. 

5. In the absence of boundary effects, for a hyperbolic PDE, neither the character of the solution, 
the error, nor the convergence rate is sensitive to the node layout. 

6. In the presence of boundaries, the solution on Cartesian nodes exhibited significant oscillations 
(Runge phenomena) near the boundary. This is not the case with hexagonal nodes, which 
was the most effective node layout in correctly capturing the physics, especially at lower 
resolutions. 

7. In all cases, quasi-uniformly scattering the nodes showed no detriment to the quality of 
the solution, error, or convergence and in the majority of the cases performed better than 
Cartesian layouts. 

8. Decreasing the viscosity by 6 orders of magnitude, (i.e. increasing in the Reynolds number by 
the same factor), does not require any change to the time step or amount of hyperviscosity 
added. 

9. In the turbulent regime, the type of node layout heavily impacts the location and structure 
of eddy development on the leading edge of the thermal as well as the degree of excitation of 
finer scale instabilities. 
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A Symmetric Stencils 


5 nodes, rad=1 9 nodes, rad=1.42 

5 I - 1 5 r 




7 nodes, rad=1 13 nodes, rad=1.74 

5 1 - 1 5 r 



-5 0 5 

37 nodes, rad=3 


43 nodes, rad=3.47 55 nodes, rad=3.61 


-5 0 5 

61 nodes, rad=4 


Figure 19: Symmetric stencils for cartesian and hexagonal nodes. Note that 13 and 37 are the only 
reasonably small stencil-sizes held in common. 
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